
clear
clear all
clear mata
clear matrix
set more off
set matsize 11000
set maxvar 30000
cap log off
capture log close
set emptycells drop
pause on


***********************************************
* USER DEFINE FILEPATH FOR REPLICATION FOLDER *
***********************************************

global path "C:\Users\wb520443\Dropbox\Research Projects\Global pollution\2_analysis\Replication"



***********************************************


global opts		a f plain coll(none) nodep nomti c(b(star fmt(%9.3f)) se(abs par fmt(%9.3f))  ) star(* .10 ** .05 *** .01) noobs nocons

global opts_fs		a f plain coll(none) nodep nomti c(b(star fmt(%9.3f)) se(abs par fmt(%9.3f))) star(* .10 ** .05 *** .01) noobs nocons


global input "$path\data"
global figdat "$path\figdat"
global figures "$path\figures"
global tables "$path\tables"




	*10 20 30 40 50
	foreach t in 10 20 30 40 50 60 70{
	    
		*local t "xx"
	    
		cap mat define countryfigure=J(1100,35,.)
			local row=1
		
		use "$input/rwi_pollution_analysis_cities.dta", replace
			drop if id==6062 & country=="MWI"
			merge 1:1 quadkey latitude longitude rwi error using "$input/ghs_vnres_rwi_table.dta", keep(1 3) nogen
			merge 1:1 quadkey latitude longitude rwi error using "$input/panel_rwi_dem_table.dta", keep(1 3) nogen
				rename std dem_std
		
			append using "$input/delhi_official_cities.dta"
			append using "$input/lagos_official_cities.dta"
			append using "$input/tbilisi_official_cities.dta"
			append using "$input/rio_official_cities.dta"

		replace error=1/error
		g both=population*error
		g both2=country_population*error
		
		g city_code2=city_code
		
		drop if city_code==.
	
		replace city_code=expanded_code
		
		replace city_code=city_code2 if city_code==.

		
		bysort city_code:egen sd_elev=sd(dem)
		bysort city_code:egen max_elev=max(dem)
		bysort city_code:egen min_elev=min(dem)
		g diff_elev=max_elev-min_elev
	
		drop if points_within_city<`t'
		*drop if points_within_city<15
		
		
		
	

		
		*drop if grid>2
		
		preserve
			keep city_code diff_elev 
			duplicates drop city_code, force
			sort city_code
			g id=_n
		
			save "$input/temp_elev.dta", replace
		restore
		
		
		levelsof city_code, local(levels) 
		
		egen total_pollution=rowtotal( power industry agri transport other)
		foreach var in  power industry agri transport other{
			g share_`var'=`var'/total_pollution
		}
		
		preserve
			keep city_code country share_* city_code2
			
			duplicates drop city_code, force
			save "$input/city_country.dta",replace
		restore
			
		//Mean, mean_ss, and max refer to how PM25 across years and months is collapse (e.g. max is the max pollution level for a RWI point across all years and months in the sample).
		foreach l of local levels {
		
				
			
			cap reghdfe mean_pm25 rwi [pweight=error] if city_code==`l', noabsorb
				cap mat def countryfigure[`row',1]=_b[rwi]
				cap mat def countryfigure[`row',2]=_se[rwi]
			cap reghdfe mean_pm25_ss rwi [pweight=error] if city_code==`l', noabsorb
				cap mat def countryfigure[`row',16]=_b[rwi]
				cap mat def countryfigure[`row',17]=_se[rwi]	
			cap reghdfe rwi dem [pweight=error]if city_code==`l', noabsorb
				cap mat def countryfigure[`row',3]=_b[dem]
				cap mat def countryfigure[`row',4]=_se[dem]
			cap mat def countryfigure[`row',5]=`l'
			sum dem if city_code==`l'
				cap mat def countryfigure[`row',6]=r(mean)
				cap mat def countryfigure[`row',7]=r(sd)
			sum rwi if city_code==`l', d
				cap mat def countryfigure[`row',8]=r(kurtosis)
				cap mat def countryfigure[`row',9]=r(skewness)
				cap mat def countryfigure[`row',13]=r(mean)
				cap mat def countryfigure[`row',14]=r(sd)
			sum mean_pm25 if city_code==`l', d
				cap mat def countryfigure[`row',10]=r(kurtosis)
				cap mat def countryfigure[`row',11]=r(skewness)	
				cap mat def countryfigure[`row',12]=r(max)-r(min)	
				cap mat def countryfigure[`row',15]=r(sd)
			sum mean_pm25_ss if city_code==`l', d
				cap mat def countryfigure[`row',18]=r(kurtosis)
				cap mat def countryfigure[`row',19]=r(skewness)	
				cap mat def countryfigure[`row',20]=r(max)-r(min)	
				cap mat def countryfigure[`row',21]=r(sd)	
				cap mat def countryfigure[`row',22]=r(mean)	
			cap reghdfe ghs_vnres_mean dem if city_code==`l', noabsorb
				cap mat def countryfigure[`row',23]=_b[dem]
				cap mat def countryfigure[`row',24]=_se[dem]	
			cap reghdfe ghs_vnres_mean dem_std if city_code==`l', noabsorb
				cap mat def countryfigure[`row',25]=_b[dem_std]
				cap mat def countryfigure[`row',26]=_se[dem_std]		
			local row=`row'+1
		
		}	
	/*	
	clear
		
		use "$input/global_pm_income.dta"
			replace ISO="IND_PAK" if ISO=="PAK"|ISO=="IND"
			collapse (mean) mean_gdp_per_cap, by(ISO)
		save "$input/temp1.dta", replace
	*/	


		
	clear
				
			svmat2 countryfigure

			drop if countryfigure1==.
			
			g country_code=_n
				run "$input/labels.do"
			
			g HI95_mean=(1.96*countryfigure2)+countryfigure1
			g LI95_mean=countryfigure1-(1.96*countryfigure2)

			g HI99_mean=(2.06*countryfigure2)+countryfigure1
			g LI99_mean=countryfigure1-(2.06*countryfigure2)
			
			
			g HI95_mean_ss=(1.96*countryfigure17)+countryfigure16
			g LI95_mean_ss=countryfigure16-(1.96*countryfigure7)

			g HI99_mean_ss=(2.06*countryfigure17)+countryfigure16
			g LI99_mean_ss=countryfigure16-(2.06*countryfigure17)
			
			foreach var in HI95_mean LI95_mean HI99_mean LI99_mean HI95_mean_ss LI95_mean_ss HI99_mean_ss LI99_mean_ss{
					replace `var'=-10 if `var'<-10
					replace `var'=10 if `var'>10
			}
			
			g tstat=countryfigure3/countryfigure4
			g tstat2=countryfigure1/countryfigure2
			g tstat3=abs(tstat2)
			
			
			
			g sig=(abs(tstat)>2.06&abs(tstat2>2.06))
			
			set scheme plotplain
				
				
			
			rename (countryfigure1 countryfigure2 countryfigure3 countryfigure4 countryfigure5 countryfigure6 countryfigure7 countryfigure8 countryfigure9 countryfigure10 countryfigure11 countryfigure12 countryfigure13 countryfigure14 countryfigure15 countryfigure16 countryfigure17 countryfigure18 countryfigure19 countryfigure20 countryfigure21 countryfigure22) (mean_beta pollution_se mean_elevation_beta elevation_se city_code average_elevation sd_elevation rwi_kurtosis rwi_skewness pm_kurtosis pm_skewness pm_range rwi_mean rwi_sd pm_sd mean_beta_ss pollution_se_ss pm_kurtosis_ss pm_skewness_ss pm_range_ss pm_sd_ss pm_mean_ss)
			
			g tstat2_ss=mean_beta_ss/pollution_se_ss
			g tstat3_ss=abs(tstat2_ss)
			
			g pm_change=mean_beta*rwi_sd
			g pm_range_share=pm_change/pm_range
			g abs_range_share=abs(pm_range_share)
			g abs_beta=abs(mean_beta)
			
			g pm_change_ss=mean_beta_ss*rwi_sd
			g pm_range_share_ss=pm_change_ss/pm_range_ss
			g abs_range_share_ss=abs(pm_range_share_ss)
			g abs_beta_ss=abs(mean_beta_ss)
			
			 
			
			
			g id=_n
			
			
			gsort mean_beta
			
			g n=_n
			
			gsort mean_elevation_beta
			
			g n2=_n
			
			gsort  mean_beta_ss
			
			g nss=_n	
			
			
	
			
			preserve
				drop if mean_beta<-5
				drop if mean_beta_ss<-5
				twoway (scatter mean_beta_ss nss if tstat3_ss>=2.06, mc(navy%80)  msize(tiny ) m(S)) ///
						(scatter mean_beta_ss nss if tstat3_ss<2.06, mc(navy%10)  msize(tiny ) m(S)) ///
							(scatter mean_beta n if tstat3>=2.06, mc(cranberry%80)  msize(tiny ) m(S)) ///
							(scatter mean_beta n if tstat3<2.06, mc(cranberry%10)  msize(tiny ) m(S) ///
							yti("Correlation between wealth and average pollution") xla(none,nogrid) xtitle("") xsc(noline) ylabel(,nogrid)  ysize(3) xsize(12) legend(off) yline(0))
							
							graph export "$figures/city_figure_`t'.png", as(png) replace
								
			restore
			
			
		
			
			
			
		
	
		
	}
		
		
	
	
	
	

	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	